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Abstract 


An analytic solution to the heavy ion transport equation in terms of 
Green s function is used to generate a highly efficient computer code for 
space applications. The efficiency of the computer code is accomplished by 
a nonperturbative technique extending Green’s function over the solution 
domain. The computer code can also be applied to accelerator boundary 
conditions to allow code validation in laboratory experiments. 


Introduction 

From the inception of the Langley Research Center heavy ion shielding program (refs. 1-3), a 
close relationship has existed between code development and laboratory experiment (ref. 3). Indeed, 
the current goal is to provide computationally efficient high charge and energy (HZE) ion transport 
codes that can be validated with laboratory integral experiments and subsequently applied to space 
engineering design. In practice, two streams of code development (laboratory codes and engineering 
c esign codes) have prevailed because of the strong energy dependence of necessary atomic/molecular 

cross sections and the near singular nature of the boundary conditions of the laboratory beam 
(refs. 4 6). 

The energy dependence of atomic/molecular cross sections is adequately dealt with by using the 
methods of Wilson and Lamkin (ref. 7) to develop efficient numerical procedures for space radiation 
(rc s. 6 and 8 10). Although these codes can conceivably be applied to laboratory validation, methods 
to control truncation and discretization errors bear little resemblance to the space radiation codes 
that need validating. Thus, a radical reorientation is required to achieve the validation goals of 
the current LaRC heavy ion shielding program, and such an approach is presented in this paper by 
extending the work of Wilson and Badavi (ref. 11). This approach consists of a nonperturbative 
technique that allows highly accurate analytic solution over the solution domain of any practical 

problem. The following section presents a brief discussion of the fundamental transport equations 
and their formal solution. 

Transport Theory 

The transport equations are derived on the basis of conservation principles. Consider a region of 
space filled with matter described by appropriate atomic and nuclear cross sections. Sketch A shows 
a small portion of the region enclosed by a sphere of radius 6. The number of particles of type j 
leaving a surface element 6 2 d(2 is given as ^ (x + 6(1, (2, E)6 2 d(2 where fc(x, (1, E) is the particle 
flux density, x is a vector to the center of the sphere, (2 is normal to the surface element, and E is 
the particle energy. The projection of the surface element through the sphere center to the opposite 
side of the sphere defines a cylinder through which pass a number of particles of type j given as 
W - 6(2, (2, E)6 d(2. This number would equal the number of particles leaving the opposite face 
if the cylinder defined by the projection were a vacuum. The two numbers of particles in fact differ 
by the gains and the losses created by atomic and nuclear collisions as follows: 


4>j(x + 6(2, (2, E)6 2 d(2 = (fffix - 6(2, (2, E)6 2 d(2 

+ 6 2 d(2 dlJ2 f cr jk (S2, (2', E, E') cf> k {x + 1(2, (2 ' , E') d(2' dE' 

— 6 2 d(2 f dl oj(E) <f>j (x + 1(2, (2, E) 
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Sketch A 

where oAE) and <r jfc (ft, ft', E,E') are the media macroscopic cross sections. The cross section 
a i (ft ft' ft, ft') represents all those processes by which type k particles moving in direc ion 
with energy ft' produce a type j particle in direction ft with energy ft. Note, several reactions can 
accomplish this result, however, the appropriate cross sections of equation (1) are the inclusive ones. 
The second term on the right side of equation (1) is the source of secondary particles integrated over 
the total volume 2 5(5 2 dft), and the third term is the loss through nuclear reaction integrated over 
the same volume. We expand the terms of each side as a Taylor series and retain terms to order 

explicitly as 

6 2 dft (x, ft, ft) + 5ft ■ V^(x, ft, ft)] = ^ d n ft ) ' m ' v< ^ (x ’ n ’ 

+ 25 W ^(ft.n'.ft.ft') d> fc (x,ft',ft') dft' dft' 

k J 

- 26a 3 (ft) 4 *] (x, ft, ft)] + 0(5 4 ) ( 2 ) 

This equation can then be divided by the cylindrical volume 25(5 2 dft) and written as 

n v ^(x,ft,ft) = E / ft', ft, ft') ^(x,ft',ft') <*ft' rfft' 

-o-j(ft) ^(x, ft, ft) + 0(5) ( 3 ) 

for which the last term 0(5) approaches zero in the limit as 5 0. Equation (3) is recognized 

as a time-independent form of the Boltzmann equation for a tenuous gas. Atomic collisions (i.e., 
collisions with atomic electrons) preserve the identity of the particle, and two terms on the right si e 
of equation (3) contribute to the gains and losses. The differential cross sections have the following 
approximate form for atomic processes: 

af k (fl, ft', ft, ft') = X] tfjn(ft') “ !) 6 Jk S ( E + en ~ E "> ^ 

n 

where the superscript at represents atomic, n labels the electronic excitation levels, and e n represents 
the corresponding excitation energies, which are small (1-100 eV m most cases) compared wit e 


particle energy E. The atomic terms can then be written as 



, n', E, E') 0*(x, ft', E') dfl' dE' - af(E) <f>j (x, n, E) 
af n (E + e n ) <f>j(x, n,E + e n ) - af{E) ^(x, ft, E) 


Tl 

« E a fn(E) 4 >M^E) + E «n E [, 7 f n (E) <j> 3 (x,n,E)] - af(E) 4 >j(x,fl,E) 

n n 

= 1e Sj ^ 


because the stopping power is 

S J (E) = '£v?n(E)en 

Tl 

and the total atomic cross section is 


(5) 

( 6 ) 


af{E) = £ af n {E) (7) 
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Equations (5) to (7) allow us to rewrite equation (3) in the usual continuous slowing down 
approximation as 


n ' V ~jE S j( E ) + a j( E ) 


4>j(x,n,E) = J ^2a jk (fl,n',E,E') 0 fc (x, «',£') dfl' dE' 


(8) 


where the cross sections of equation (8) now contain only the nuclear contributions. For convenience, 
we can rewrite this equation as the following matrix equation: 


B0 = P0 (9) 

where B is a diagonal differential matrix operator consisting of the drift term, energy loss term, 
and collisional loss term on the left side of equation (8); 0 is a vector of the particle fields with the 
heaviest particles in the lowest entries; and P is the collisional source matrix that is assumed to be 
an upper triangular integral operator. 

Equation (9) may be rewritten as 


0 = 0 O + B ! P0 (10) 

where 0 O is the solution of the homogeneous equation and is equal to the incoming flux at the 
boundary. A propagation matrix Green’s function G 0 , which is the inverse of the multivariate 
integrating factor (ref. 3) for the right side of equation (9), is defined as follows: 

0 = G O 0 B + B _1 P0 (11) 

where 0# is the boundary condition and 


BG 0 = 0 (12) 

The Neumann series solution of equation (11) can then be written as 

0 = G 0 0 B + B -1 P G 0 0 B + B -1 P B -1 P G 0 0 B + . . . (13) 


3 



The series converges after N terms where N is the rank of the space because B 1 P G 0 is a contraction 
mapping operator which is upper triangular (ref. 2). The complete Green s function is defined as 

0 = G0 B (14) 


and is given by the series 


G = G„ + B -1 P G 0 + B _1 P B l P G 0 + ... 


(15) 


and satisfies 

G = G 0 + B _1 PG (16) 

where B 'P is factored from all but the first term on the right side of equation (15). Clearly, G 
satisfies the same operation equation as 0 in equation (9). Thus, 

BG = PG (17) 


where G is a matrix operator that reduces to the identity at all points on the boundary. In component 
form, equation (17) is written as follows: 


n ■ V - Tfg Sj( E ) + Vj G jm(E 


,Eo, n,x) = W a jk (E,E',n,n') G km {E\E 0 ,n\x) dtf dE ’ 

(18) 


where Sj(E ) is the change in E per unit distance, and E now denotes the energy per nucleon. An 
arbitrary solution to the Boltzmann equation within a closed convex region can be written as 


^•(£,n,x) = z / Gjm(£,E',n-n',x-r) * f m (E’,n\ r) drt dE' dr 


(19) 


where fm{E f } il\T) is the incident flux at the boundary (ref. 12). Because the transport problem is 
formulated in terms of a single Green’s function algorithm, the validation of the Green s function in 
the laboratory meets our objective of having a space- validated code (refs. 13 and 14). 

The first step is to develop an equivalent Green’s function algorithm in one dimension that 
matches our current capability in space radiation shielding (refs. 6 and 15). The algorithm is based 
on the closed-form solution to the following one- dimensional equation for a monoenergetic beam at 
the boundary: 

(l-A §,(£) + <3 lm (E,£„.x) = £/ a jk (E,Ef) G km tf,E 0 ,x) dE' (20) 


The possibility of validation for 20 Ne beams of this algorithm (with multiple scattering corrections) 
has already shown promise (refs. 5 and 16), but improvements in the nuclear data base are required 
to achieve agreement. If we restrict our considerations to multiple charged ions, then the right side 
of equation (20) can be further reduced to 


d_ 

dx 




Gjm{E,E 0i %) — } ^ &jk E 0 , x) 

k 


( 21 ) 


for which we now consider the solution. 
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Approximate Green’s Function 

We now derive an approximate Green’s function from equation (21), which can be simplified by 
transforming the energy into the residual range as 


r E d&_ 

'o Sj(E') 


and defining new field variables as 


if>j(x,rj) = Sj(E) <f>j(x,E) 


( 22 ) 


(23) 


Gjm{x, rj,r' m ) = Sj(E) G jm (x, E, E ') (24) 

so that equation (21) becomes 

Qj. 4" Oj J ~ ^ 1 (26) 

where i/j is the range scaling parameter, Uj = Zj/Aj. The boundary condition for the Green’s 
function is 

Gjm(Qi r jt r m) ~ fijm ^{ r j ~ r m) (26) 

and the full solution for an arbitrary boundary condition is 

roc 

~ ^ ^ I Gjm{ x i ^m) 

m ^ 


The solution to equation (25) can be written as a perturbation series as 

Gjmi x y r j > r m) = r j’ r m) 


where 


and 


r j> r 'm) = 9(j) Sjm S(x + Tj - r' m ) 


0jm( x ’ r > r m) 


Vj &jm 9(ji 




x{y m ~ vj ) 


where Q Tj, r^) is zero unless 


■'jm 


1/j , x , Vj 

T7~( r j + x) < r m < — r j + x 


The second-order term is 


and is nonzero for 




Vjk a km 9(3 , k, m) 

r > _ r ; 

r ju r jt 


r ml — r m — r m 


(28) 

(29) 

(30) 

( 31 ) 

(32) 

(33) 
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where 


' v i 

± r 3 +X 

(Vm > V k - > Vj) 

v jTj + u k x 
Vm 

{v k >Vm> Vj) 

. Sr/j + x 

(l/ w > Vj > V k ) 

% ( r J + x ) 

{v m > V k > Vj) 

£<»■>+*> 

(l/ k > I'm > Vj) 

VjTj + Vk X 
Vm 

( v m > Vj > v k ) 


and Tj u and r f j £ are found at the limits of equations (34) and (35). The third-order term is 

^,(3)/ / \ a jk a ki ®tm > m ) /o£\ 

r V r m) » Zv T~Z7' ( 36 ) 

k } £ 3 U J? 

Higher order terms are similarly derived. In equation (36), the function g of n arguments are given 
by 

g(j) = exp (-<Tj x) (37) 


9(j J '2 j * • * * 3ni in+l) ““ 


9(j 1 ? J 2 ? Jn-liJn) • - - > 3n-h jn+l) 


a 3n+\ a 3n 


As noted in reference 11, the solution to equation (21) can be written in terms of equation (38) as 
follows: 

ipj(x,rj) = exp {—cr j x)fj(rj + x) + ]T G^{x) ( F m {r' mu ) - F m (r^)j (39) 


(i) . . <7 j 3 i <7 JiJ 2 ’ ” * » a jn- 2 m 9{3i3i32> • - ' i Jn- 2, m ) 

- 2s A (i) 

3 1 '3 2 ’ *■*> 3 n -2 


a(D=x 1^-1 


for (z = 1) 


and for (z > 1) 

x _ 0 > Uk > 

A w = < z - l) (v k > v m > Vj) > 

, x fe - ?) (l/m > > ^ > 

In equation (39), F m (r) is the integral flux at the boundary and can be represented as 


/ CX) 

fm(r) dr ' 


Implementation of equation (39) can now be accomplished independent of the character of the 
boundary values / m (r ; ) to give accurate results for both space and laboratory applications. 
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Results of Perturbation Series 


The solutions given in equations (23) to (43) contain a finite number of terms, but higher terms 
contain large numbers of components. Thus, only a few terms are usually generated. To better 
understand the nature of the solutions given by equations (30), (32), and (36), we now examine 
the perturbation series approximation to the scaled Green’s function for a broad set of fragments. 
Figure 1 shows the first three collisional contributions to the Green’s function for producing 52 Cr. 
In the figure, the atomic weight of each projectile (A p ) is truncated to the nearest isotope. The main 
contribution comes from near projectile elements up to third collision terms. Figure 2 shows the 
relative contribution of each term. The third collision term contributes less than 20 percent over the 
domain, and the contributions of higher order terms are smaller. Although fourth collision terms 
are not negligible at 50 cm, the total 52 Cr production at this depth is small. (See fig. 3.) 

Figures 4 to 6 show similar results for 40 Ar production. The importance of the higher order terms 
in producing the lighter fragment is apparent when figure 2 is compared with figure 5. For example, 
the third collision for 40 Ar is as large as 40 percent in the displayed domain. Figures 7 to 9 show the 
28 Si production terms. The third-order 28 Si production term dominates at large depths, and higher 
order terms are expected to be large. The main 28 Si production comes from collision with primary 
ions up to mass 42. Figures 10 to 12 show similar results for 16 0 production. Clearly, terms of fourth 
and higher orders must be added into any meaningful calculation, but such terms require lengthy 
computational procedures. To evaluate the Green’s function to all orders of collisions, we now derive 
a nonperturbative procedure. Thus, convergence problems associated with the calculations described 
in this section are circumvented. 


Nonperturbative Procedures 


In generating the graphs in figures 1 to 12, we found that the large number of terms contributing 
at large depths leads to an inefficient computational procedure. One of the difficulties is that many 
isotopes can be produced in successive collisions of the incident ions and their collision products. 
The first three terms of equation (28) are evaluated rather rapidly; thus, one question remains. How 
are the higher order terms to be made efficient? The following example is one possible approach. 

The generator equation of the ^-functions is given as a solution to the following equation: 


d 1 


9jm( X ) = E a 3k 9 krn (x) 


(44) 


This equation is subject to the boundary condition gj m ( 0) = 6j m . Thus, it can be inverted as 

J rx 

l E a jk exp [-ajix - y )] g k Jy) dy (45) 

0 k 


where m denotes the solution for incident ion type m at the boundary (x = 0). Assume the solution 
of this equation is 

gfrl(x) = &jm exp {-(J 3 x) (46) 

and the iterated equation is 



E a jk exp [— <7j(x — y)] gi m l) (y) dy 

k 


(47) 
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The higher order terms can be summed as 


E f exp \—(T J {x — y)] 

i>2 J ° k 


9kl(y ) + E &y) 

i> 2 


dy 


= 9?m(x)+ [ E ° ik exp [—Oj{x - y)] E sSniv) d V 

Jo k i > 2 


( 48 ) 


We find the following relation for the full series: 

9 jm ( x ) = E 9 ]k ( x - y) 9 k Jy) 

k 


(49) 


The appendix contains a derivative of this equation. Note that just as g jm (x) is the solution of 
equation (44) with g jm ( 0) = 6 jm , 9 jm ( x ) is the solution of equation (44) having the value g jm {y) at 
x = y. Equation (49) can be used to generate solutions to all orders and arbitrary depths. Because 
equation (49) is an algebraic relationship that is true for an arbitrary y, we expect a rapid evaluation 
to arbitrary order and highly efficient computations, especially on vector processors. It follows from 
equations (29), (30), (32), (36), and (40) that the full Green’s function can be written as 


9jm ( x ) — 9(3) bjm a jm m ) 


A( 2 ) 


g jm (x, rj,r' m ) = exp (-(Tjx) S jm 6{x + rj - r' m ) + 9(j > m ) + _ 

(50) 

where the domains of successive terms are given by equations (31) and (33). Equation (50) can now 
be used to generate solutions to the ion transport problem. 


A practical computation consists of three steps. The first step is to use the recurrence relation in 
equation (49) to adequately cover the solution domain with g jm (x). The second step is to evaluate 
the differential and integral spectra at the boundary by equation (43), and the third step is to use 
the following equation to evaluate the solution over the domain: 


ipj(x, rj) = exp (-ajx)fj(rj + x) + E ~ l Fm ^ r mu) ^m(r m ()] 

+ £ (F m (C„) - F„( r^>] 

m,k ^ 


(51) 


where r' mu and r' m( in the second term are given by the inequality limits of equation (33), and r mu 
and r' m( in the third term are given by equations (34) and (35). The quantity denoted as 9 jkm (x) is 


given as 

9 ]km ( x ) = 9 Jk ( x - y)9 km (y) 


(52) 


These nonperturbative techniques hold great promise for accurate and efficient computational 
methods for evaluation of the HZE particle fields in space and laboratory problems. However, the 
techniques must still be extended to light ion and especially neutron fields. 


Results of Nonperturbative Procedures 

A computer program has been written to calculate the Green’s function with these nonpertur- 
bative methods. Figure 13 shows the values for the collision related terms of Gjm{ x ^ r j,r' m ). In the 
figure, x is the depth in a water medium, and Z p is the charge of the incident projectile. Specif- 
ically produced species are noted in the figure. The production of any given species is dominated 
by the projectiles of nearly the same but greater charge. Figure 14 shows the relative contribution 
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of the first collision term. The multiple collision terms are more important for those projectiles 
whose charge is far removed from the specific species. For example, the first collision term con- 
tributes ^20 percent to the flux at 50 cm for Ni ions ( Ap = 59) to a third-order approximation 
(fig. 11(a)), however, it only contributes about 10 percent to the O flux when all orders are consid- 
ered (fig. 14). The Green s function and equation (51) are used to evaluate the composition of a 
monoenergetic 535 MeV/nucleon Fe beam in a water column at several depths; figure 15 shows the 
results. These results will be compared with other results from experiments in the near future. The 
greatest advantage of the nonperturbative method is the computation efficiency. For example, to 
complete the computation for figures 1 to 12 required a week of time on a Digital Equipment Corp. 
VAX 785 compared with a few minutes to complete the computations for figures 13 to 15. 

Concluding Remarks 

The Green s function is a consistent method for developing space shielding codes that can be 
validated in laboratory experiments. Perturbation theory is a means of generating the Green’s 
function for such purposes, but it is hampered by the inefficiency of the computational procedures. 
A nonperturbative technique is derived in which high computational efficiency is achieved without 
loss of accuracy. These nonperturbative methods hold great promise in allowing highly efficient 
computer codes for space shielding. The computer code can also be applied to accelerator boundary 
conditions to allow code validation in laboratory experiments. 


NASA Langley Research Center 
Hampton, VA 23681-0001 
January 22, 1993 
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Appendix 

A Derivation of Equation (49) 

The generator equation (44) for g jm (x) has the same form and boundary conditions as the Green’s 
function equation (17). The corresponding particle fluence equation (9) is given by 

(t£r + CT? ) = '^2 a J k ^ k ^ 

where <PkW) ' H known, and the Green’s function solution (14) is given by 

*,(i) = E »*.(*> ^<°> (A2: 


We also have 


My) - ^2 9km(y) MW) 

m 


where y denotes a point on the J^-axis. 
For x > y, we can also write 


4j(x) = ]C 9 jk ( ' X - ^ fc ( y ) 

k 


Substituting equation (A3) into equation (A4) gives 

4>j(x) = 3jk( x ~ y) 3k m (y) MW) 

m k 

By comparing this equation with equation (A2), we obtain the result 

g jm (x) = J2 - y) 9 ^(y) 


(A3) 


(A4) 


(A5) 


which is equation (49). 
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(a) First collision. 


(b) Second collision. 


(c) Third collision. 


Figure 4. First three collision contributions to 40 Ar production. 



(a) 


First collision. (b) Second collision. (c) Third collision. 

Figure 5. Fractional contribution to 40 Ar production from first three collision terms. 



Figure 6. Total 40 Ar production from first three collisions. 
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